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Abstract 



In the diffusion-collision model, the unfolding or backward rates 
are given by the likelihood of secondary structural cluster dissociation. 
In this work, we introduce a backward rate calculation modeled from 
a Kramers-type thermal tunneling through a barrier, which represents 
the free energy potential well for buried hydrophobic residues. Our 
results are in good agreement with currently accepted values and the 
approach suggests a link between the diffusion-collision and folding 
q funnel models of protein folding. 



I. Introduction 

In the diffusion-collision model (DCM) of protein folding (Bashford et al 
1988, Karplus and Weaver, 1976, 1979, 1994) the protein is modeled using 
a collection of spheres connected by flexible strings. The spheres represent 
the secondary structural elements such as a-helices or /3-sheets, called mi- 
crodomains, that make up the protein. 
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The folding process from a completely unfolded protein to the the final native 
state is accomplished via diffusion through the solvent, collision, and finally 
coalescence of the microdomains. The state of the protein is defined by the 
number of pairings (hydrophobic interactions) between the microdomains 
that are present at a given time t. The rate equations for transitions between 
these states can be written as 

where P(t) is the vector of states and K is a matrix containing the transition 
rates between the different states. 

As an example, consider a simple monomeric protein containing two a-helices 
A and B joined by a connecting string. This gives us a one-pair/two-state 
system. We will call state 1 the state when the two microdomains are not 
in hydrophobic contact, and state 2 when they are hydrophobically docked. 
In this case (||) can easily be written out explicitly 



A(t) = -h^ 2 Pi(t) + k 2 ^ 1 P 2 {t) 

h(t) = h^Pxit) - k 2 ^iP 2 (t) 

or in matrix form as 

d_ ( Px(t) \ 
dt [ P 2 (t) J 

A more complicated protein having, say, n microdomains would involve 
p = n{n — l)/2 pairings, 2 P states Pi(t) and a 2 P x 2 P rate matrix K. In 
general the calculation of the elements of the rate matrix K is somewhat 
involved. The forward rates are the rates of microdomain coalescence. In 
the diffusion-collision model the forward rates are calculated assuming the 
microdomains diffuse through a solvent environment, the space of which is 
limited by the length of the intervening strings and the van der Waals radii 
of the microdomains. These microdomains are assumed to be nascently 
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formed, and their degree of formation is given by a helix-coil transition the- 
ory (Zimm and Bragg, 1959, Lifson and Roig, 1961) as in AGADIR (Lacroix 
et al, 1998, Munoz and Serrano, 1994 a,b,c, 1997) calculation in the case of 
a-helices, or via a combination of theory (Munoz et al, 1998) and experiment 
(Dinner et al, 1999) in the case of /3-sheets. As the microdomains undergo 
diffusion, they occasionally collide. When this happens the microdomains 
coalesce with a probability 7, being held together by hydrophobic interac- 
tions in the case of a-helices, or a combination of hydrophobic and hydrogen 
bond interactions in the case of /3-sheets. The coalescence probability 7 is 
given by the likelihood that the microdomain is in a-helical or /3-sheet form, 
the percentage of hydrophobic area, and the likelihood of proper geometrical 
orientation upon collision. 

The forward folding times in the mean first passage time approximation 
(Weiss, 1967, Weaver, 1979, Szabo et al 1980) are given by 

f D 7 DA { ' 

where V is the diffusion volume available to the microdomain pair, A is the 
target area for collisions, D is the relative diffusion coefficient, 7 is the prob- 
ability of coalescence upon collision and I and L are geometrical parameters 
calculated for diffusion in a spherical space. The inverse of the first pas- 
sage time-scales r/ are the forward folding rates kf that are used in the rate 
matrix K . In the example given by (§) and (||), kf is &i_>2- 
The microdomain pairings can also dissociate. To date, in typical diffusion- 
collision model calculations, the form of the backward folding, or unfolding 
times 17 used for two microdomains A and B comes from from the Van't 
Hoff-Arrhenius law (Van't Hoff, 1884, Arrhenius, 1889) given by 

Tb = y~\^ABlk B T (5) 

where / is the free energy change per unit buried hydrophobic area in the 
pairing (Chothia, 1974), Aab is the buried area (Lee and Richards, 1971), 
ks is Boltzmann's constant, T is the temperature and v is an attempt rate. 
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In diffusion collision model calculations to date the attempt rate is a param- 
eter whose value must be set by hand, usually requiring some guess work or 
adjustment to obtain the desired result (see for example Burton et. al 1998). 
The typical values used lie in the range 1 x 10 9 s -1 — 1000 x 10 9 s _1 , obtained 
from estimates of covalent bond oscillation frequencies (Fersht, 1999). This 
procedure is not explicitly stated in most diffusion-collision model studies; 
it is justified, however, when the equilibrium occupation probabilities of the 
states are known. In fact, the rate constants can be calculated from such 
probabilities (Chandler, 1987). In cases where the final occupation proba- 
bilities are unknown, for instance in studies of protein mis-folding and non- 
native kinetic intermediates (Beck et al., 2000), such methods are clearly not 
possible. 

In this work we show how the parameter v, or more generally the unfold- 
ing rates, can be determined from thermal fluctuations providing a means 
of avoiding the guesswork. This makes the diffusion-collision model more 
predictive and enables it to be used in situations where the final occupation 
probabilities are unknown. 

II. Calculation of the Unfolding Rates 

In order to calculate the backward rate it is convenient to view the unfold- 
ing process as that of diffusion within, and escape from, an effective one- 
dimensional potential well. This is a good approximation if only motion per- 
pendicular to the hydrophobic contact surface is important in microdomain 
pair dissociation. 

Consider the pair of microdomains A and B connected by a string. Mi- 
crodomain coalescence and dissociation can be approximated by diffusion in 
a potential like the one depicted in Figure 1. The quantity that diffuses is 
the separation between the microdomains and we can think of 

dP( x ) = p{x)dx (6) 
as the fraction of microdomains with a distance between them in the range 
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STATE 2 STATE 1 



Figure 1: Diffusion potential for the two microdomains A and B. The poten- 
tial is infinite on the far right because the microdomain separation cannot 
exceed the length of their connecting string. On the left it is also infinite, 
in this case because of the hard-core repulsion of the van der Walls con- 
tact between the microdomains. The potential energy barrier has a height 
Eb = /Aab, the free energy difference between paired and unpaired states 
with a buried hydrophobic area Aab, and width L, which we take to be the 
diameter of a water molecule. 



x and x + dx. In a potential V(x) the probability density p(x) satisfies the 
one dimensional Smoluchowski equation 



dp(x,t) d 
dt dx 



dx 



dx 



(7) 



The outer boundary on the far right is given by an infinite potential, arising 
from the two microdomains coming to the end of their connecting string. 
The inner boundary is also infinite, stemming from the hard-core repulsion 
of the van der Waals contact between the microdomains. The well depth 
Eb = f Aab is the free energy difference between paired and unpaired mi- 
crodomains and the well width L is set to the diameter of a water molecule. 
A separation larger than L exposes the microdomain to the solvent, the free 



5 



energy savings is lost, and the microdomains separate, resulting in an escape 
from the potential well. 

The forward folding rates can be satisfactorily found by the mean first pas- 
sage time using the method outlined above (|). As a consequence of this we 
shall forego direct analysis of the entire diffusion process and concentrate on 
the escape of a bound microdomain out of the potential well in order to find 
the backward folding rate. 

If we concentrate on the part of the potential where x < L in Figure 1, we can 
think of the potential as being perfectly reflecting on one side and partially 
permeable on the other. The permeability arises from the thermal fluctu- 
ations in the energy which allow the microdomains to escape the bounded 
region, provided they have sufficient energy, and it can be expressed using 
an appropriate choice of boundary conditions. We proceed as follows. 

For the inner boundary, representing van der Waals contact, we impose per- 
fectly reflecting boundary conditions, meaning the flux through the boundary 
is zero 

ax 

The flux at the permeable boundary at x = L depends on the density of 
microdomains at that boundary and the probability that the microdomain 
energy is high enough to thermally tunnel through the boundary. The dif- 
ferential element of flux at the boundary L of microdomains with relative 
velocity between v and v + dv is given by 

dJ(L,i) =vp(L,t)dN(v) (9) 

where 

dN ^ = {dsi) l ' 2e ^ nkBTdv (10) 

is the Maxwell-Boltzmann distribution, namely the fraction of microdomains 
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with relative velocities between v and v + dv and fx is the reduced mass given 
by 

m A m B 

A* = 7 \ (11) 

(771^+777,73) 

where tua and tob are the masses of the two microdomains. The distribution 
is normalised so that the bounds of the integral go from — oo to oo. 
It is clear that in order to find the total flux through the outer boundary at 
L we must integrate over all microdomain velocities larger than + >/ Eb/2m 
since the potential well height is Ej,, and only microdomains with velocity 
higher than that can escape and thus contribute to the flux leaving the well. 
This yields a flux out of state 2 

jr\L,t) = p(L,t) (^) 1/2 ^ bAsT (12) 

If we assume that the probability is uniformly distributed across the well, 
then p(x,t) = P2(t)/L everywhere within the well. This is a reasonable 
assumption because the permeability of the well is so low. Furthermore it 
can be shown that typically the damping time (or velocity autocorrelation) 
time is considerably smaller than the well crossing time and the mean first 
passage time. Clearly then, diffusion in the interior of the well distributes 
the microdomain separations evenly among the available space and therefore 

j r(L,t) = ^^) 1/2 e- E ^ T . (13) 

The interesting result is that by using a stationary flux method (Kramers, 
1940), or equivalently setting = 1> the outward flux at the boundary 

in one dimension J out (L,t) is the fraction of microdomain pairs crossing the 
boundary at L per unit time and thus clearly identical to the unfolding rate 
kb for the AB microdomain pairing. Below we explain this in more detail 
and show more generally how the rate equations of the diffusion-collision 
model ([j], |2|) can be derived from the diffusion equation. 
To understand this, we remember that the Smoluchowski equation ([?]) is 
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simply a detailed statement of the equation of continuity 



^ = V.J(x, t ) (14) 

The quantity in brackets in ([?]) is the flux J(x, t), which includes the ordinary 
Vp(x, t) term, and the second term which takes into account external forces. 
Integration of both sides of over the bound part of diffusion volume, 
using the divergence theorem on the R.H.S. and assuming isotropic flow into 
and out of the boundary of the diffusion volume yields 



dt 



n-J(t)A = {Jl n (t)-J^(t)}A. (15) 



The LHS of this equation is the rate of change of probability in the bound 
region. On the RHS n is a unit vector normal to the boundary, n • 3(t) is 
the net flux which traverses the boundary surface and A is the area of that 
boundary surface. We suggestively write the net flux n • J(i) in terms of 
J 2 n (t) and J 2 (£), the probability fluxes in and out of the diffusion volume 
in question. 

We see from comparison to @ that (15) can be written out explicitly for 



our case 



{Jf(i) - JT\t)}A = k 1 ^ 2 P 1 (t) - k 2 ^P 2 (t). (16) 



One can easily see that the positive terms on either side of (|16|) are equal to 
each other as are the negative terms. Referring to (111), this quantity is 



r 2 u \t)A = k 2 ^p 2 {t) (17) 

making the rate 
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(18) 



In this expression the one-dimensional "area" through which the probability 
is flowing is simply a point and therefore A = 1 so identifying Jlp^ty) with 
(|l3|) we can solve for the unfolding rates in K, &2->i in (H) and (§), by dividing 
( |i~3| ) by P2{t), or setting P2(t) = 1 in ( filf ) according to the stationary flux 
method. We find the backward folding rate in the one-dimensional case to 
be 

The terms preceding the exponential correspond to our prediction for the 
Van't Hoff-Arrhenius attempt rate v in (||). As an example, the attempt 
rate found for a coalesced pair of 16-residue Regan-Degrado (Regan and 
Degrado 1988) helices with a combined hydrophobic area loss of 600A 2 is 
64 x lO 9 * -1 . 

We believe that the most probable dissociation of a microdomain-microdomain 
pairing occurs via relative motion along a vector connecting the centers-of- 
mass of the two microdomains. This implies that typically the one dimen- 
sional case is the best way to view the dissociation event and therefore that 
one should use (19) to go about calculating the rate. The one dimensional 



case should also be sufficient if there is an "unzipping" of paired a-helices. 

It is possible, however, that dissociation may also include a relative rolling 
motion or other motion perpendicular to the axis connecting the microdomain 
pairs. In this case one needs to repeat the calculation above with a few mi- 
nor differences: The relative velocity distribution of the microdomains is still 
the one dimensional Maxwell-Boltzmann distribution because the degrees of 
freedom parallel to the surface through which the probability is flowing do 
not contribute to escape from the well, the probability in the bound region 
is assumed to be evenly distributed in a two-dimensional "volume", namely 
p(x,t) = P2(t)/iTL 2 in the two-dimensional analogue of (|T^), and the flux 
goes through a two-dimensional "area" A = 2ttL in (111). This calculation 
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yields the result 



»" = !(^) e_w ' (20) 

Due to the steric clashing of the side chains it seems rather unlikely that 
dissociation would include a relative sliding motion along the axes of the 
microdomains. For completeness, however, we include the three-dimensional 
result derivable from analogous considerations to the ones given above 



t2 - 1 = I l2^J ' " ' ■ (21) 

Although the three-dimensional result seems an unlikely candidate for pro- 
tein unfolding it may be relevant in the context of molten globules. 

This approach succeeds in removing the free parameter v from the model, 
and allows us to find the backward rates from a simple energetic model based 
on diffusion in a potential with appropriate boundary conditions. Moreover 
our results for the one, two and three-dimensional unfolding rates have a 
\JT ' 1 1 fx dependence that could be used to distinguish between this and other 
proposals for the mechanism of helix-helix dissociation. 

The removal of the parameter v is important when considering folding pro- 
cesses which do not involve the native state. In previous applications of the 
diffusion-collision model, the folding kinetics from a denatured or random 
coil state to the final native state were followed. In such a case, it is rea- 
sonable to set the parameter v such that the native state achieves 90 or 95 
percent of the probability, because we know that the final state is attained 
at the end of the folding process. In studying intermediate processes or more 
importantly, non-native intermediates (Beck et al, 2000), where the occupa- 
tion probability may be completely unknown, such reasonable estimates of v 
are not available. In this case, elimination of v as a free parameter is crucial. 
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III. Concluding Remarks 



We have presented a calculation for the helix-helix dissociation rate using a 
simplified potential surface, a square potential having a depth equal to the 
free energy savings of hydrophobic docking and width equal to the diameter 
of a water molecule. We have found the unfolding rates arising from thermal 
fluctuations out of this potential well to be in good agreement with currently 
accepted values of v. The potential itself is due to hydrophobic forces, which 
to date are not well understood and for this reason the potential chosen was 
simple. 

The initial motivation of this work was to eliminate the free parameter v 
from the diffusion-collision model. In the context of this model our result 
should allow us to perform kinetics simulations than were not possible before, 
for instance time developments of non-native intermediates, in which the 
occupation probabilities are unknown, and reasonable estimates for v are 
not possible. 

The results presented here also predict a u oc \JT j [i dependence in all cases 
which can be distinguished experimentally from other proposals such as the 
covalent bond model where v oc ksT/fi (Fersht, 1999). Another difference is 
the dependence of the unfolding rates on the states, not only through the hy- 
drophobic area, but also through the reduced mass \x of the microdomains or 
groups of microdomains undergoing dissociation. This is markedly different 
from typical diffusion-collision model calculations where the attempt rate v 
is assumed to be the same for all dissociation events within the protein. 

Along the way we have re-examined the idea that couching the problem 
of association and dissociation of microdomain pairs via diffusion over a 
potential surface affords us a clear and simple picture of the protein folding 
process. Indeed, the rate equations of the diffusion-collision model can be 
derived from such a picture. Interestingly, the potential energy surface is 
an element of other models such as folding funnels (for an overview see 
Bryngelson et al 1995). 

Admittedly, our approximation is rather rough. The thrust of future work 
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should therefore be in constructing more realistic potentials including acti- 
vation energies of hydrophobic dockings. 

Generalization to more complex proteins is straightforward because every 
interaction between helices or clusters of helices can be considered a two- 
state process similar to the one described above; generally, the folding or 
unfolding of a given protein involves several such processes. In this case 
the diffusion space for each possible pairing could be constructed, and the 
forward and backward rates for each transition can be found, as outlined 
above, to construct the transition rate matrix K in order to find the time 
evolution of the state vector P(i). This is a much more tenable proposition 
than directly solving the Smoluchowski equation (|7|) on such a complicated 
potential surface, although it could be done in principle. 
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